% Thus function solves the DSGE model under perfect foresight by solving a
% fixed point problem using the standard Newton-Raphson routine.
function [Z,f,diff,iter,errorMes] = FixedPointSolver(Z0,x_t,y_tN,N,setupEPer)

% Unfold setupEPer
ny       = setupEPer.ny;
mx       = setupEPer.mx;
myx      = setupEPer.myx;
dispOn   = setupEPer.dispOn;

diff     = 1;
diffOld  = 100;
Zvec     = reshape(Z0,N*(ny+mx),1);
f        = DSGEforesight_N(Z0,x_t,y_tN,N,setupEPer);
iter     = 0;
errorMes = 0;
updateJ  = 1;
stopflag = 0;
invW     = zeros(ny+mx,ny+mx,N);
fx11     = zeros(ny+mx,mx,N);
fx12     = zeros(ny+mx,myx,N);
fyp      = zeros(ny+mx,ny,N);

while abs(diff) > setupEPer.tolf && iter < setupEPer.MaxIter && stopflag < 2
    iter = iter + 1;
    % The gradient
    if updateJ == 1
        if setupEPer.JacobianOption == 1
            epsValue = 1e-8;
            dimf     = N*(ny+mx);
            fp       = zeros(dimf,dimf);
            fm       = zeros(dimf,dimf);
            for i=1:dimf
                % Positive step
                ZvecPlus      = Zvec;
                ZvecPlus(i,1) = Zvec(i,1) + epsValue;
                fp(:,i)       = DSGEforesight_N(reshape(ZvecPlus,ny+mx,N),x_t,y_tN,N,setupEPer);
                % Negative step
                ZvecMinus      = Zvec;
                ZvecMinus(i,1) = Zvec(i,1) - epsValue;
                fm(:,i)        = DSGEforesight_N(reshape(ZvecMinus,ny+mx,N),x_t,y_tN,N,setupEPer);
            end
            % Gradient
            J = (fp - fm)/(2*epsValue);
            ZvecNew  = J\(-f) + Zvec;
            % For debugging
            %Janal = analJacobian(reshape(Zvec,ny+mx,N),x_t,y_tN,N,setupEPer);;
            %Janal = full(Janal);
            %DIFF = round(Janal-J,6);
            %max(max(abs(Janal-J)))
        elseif setupEPer.JacobianOption == 2
            J = analJacobian(reshape(Zvec,ny+mx,N),x_t,y_tN,N,setupEPer);
            ZvecNew  = J\(-f) + Zvec;            
        elseif setupEPer.JacobianOption == 3
            [delta,invW,fx11,fx12,fyp] = nextStepFixedPointSolverBoucekkine(reshape(Zvec,ny+mx,N),x_t,y_tN,N,...
                reshape(f,ny+mx,N),invW,fx11,fx12,fyp,updateJ,setupEPer);           
            ZvecNew = delta + Zvec;
        end
    else
        if JacobianOption == 1 || JacobianOption == 2 
            ZvecNew  = J\(-f) + Zvec;
        elseif JacobianOption == 3
            [ZvecNew,invW,fx11,fx12,fyp] = nextStepFixedPointSolverBoucekkine(reshape(Zvec,ny+mx,N),x_t,y_tN,N,...
                reshape(f,ny+mx,N),invW,fx11,fx12,fyp,updateJ,setupEPer);  
        end
    end
    if dispOn == 1
        tmp = reshape(ZvecNew,ny+mx,N);
        subplot(2,1,1)
        plot(tmp(1:ny,:)');
        title('controls')
        subplot(2,1,2)
        plot(tmp(ny+1:ny+mx,:)');
        title('states')
        figure(gcf);
    end
    diff = max(abs(ZvecNew - Zvec));
    
    % Updating
    if diff < diffOld 
        f        = DSGEforesight_N(reshape(ZvecNew,ny+mx,N),x_t,y_tN,N,setupEPer);
        Zvec     = ZvecNew;
        stopflag = 0;
        diffOld  = diff;
        if setupEPer.updateJDefaultOn == 1
            updateJ  = 1;
        else 
            updateJ = 0;
        end
    else 
        updateJ  = 1;
        stopflag = stopflag + 1;
    end
    % We stop if the residuals of f0 are smaller than residualMax
    if setupEPer.residualMax > 0
        if max(abs(f(1:ny+mx,1))) < setupEPer.residualMax
            diff = 0;
        end
    end
    if dispOn == 1
        disp(['Change in Z = ', num2str(diff)])
        disp(['Udate if J  = ', num2str(updateJ)]);
        disp(['Stopflag    = ', num2str(stopflag)]);
        disp(' ');
    end
end
if setupEPer.residualMax > 0
    % For the hybrid case where the criteria implied by setupEPer.residualMax is sufficient.
    if abs(diff) > setupEPer.tolf || iter > setupEPer.MaxIter || any(isnan(diff))
        errorMes = 1;
    end
else
    % Not for the hybrid case
    if abs(diff) > setupEPer.tolf || iter > setupEPer.MaxIter ||  any(isnan(diff)) || max(abs(f)) > setupEPer.tolf
        errorMes = 1;
    end
end
Z = reshape(Zvec,ny+mx,N);
end

